Thermal roughening of an SOS-model with elastic interaction 
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We analyze the effects of a long-ranged step-step interaction on thermal roughening within the 
framework of a solid-on-solid model of a crystal surface by means of Monte Carlo simulation. A 
repulsive step-step interaction is modeled by elastic dipoles located on sites adjacent to the steps. 
In order to reduce the computational effort involved in calculating interaction energy based on 
long-ranged potentials, we employ a multi-grid scheme. As a result of the long-range character of 
the step interaction, the roughening temperature increases drastically compared to a system with 
short-range cutoff as a consequence of anti-correlations between surface defects. 



I. INTRODUCTION 

At low temperatures crystal surfaces are known to as- 
sume the shape of a plane facet. With increasing temper- 
ature fluctuations gradually contribute a nonzero thick- 
ness to the initially flat facet. This surface thickness 
finally diverges at a finite temperature, the roughen- 
ing temperature, where the order of the facet is lost 
completely. This transition can be described by a set 
of renormalization group equations first analyzed by 
Kosterlitz and Thoulesil. Because of its unusual prop- 
erties and the relation to the two-dimensional Coulomb 
gaso, thisjjOjughening transition has attracted substantial 
attentionfflit 

Various discrete solid-on-solid models have been shown 
to undergo this type of transition. Most of these mod- 
els incorporate local interactions, at most next-nearest 
neighbor interactions. Within some of these models a 
transition involving in-planc dkp|ipd^is-|papible, usually 
referred to as prerougheningflffllMjyy'By. Step-step 
interaction by means of elastic deformation of the crystal, 
however, is of a long-ranged nature and has apparently 
not been previously studied in the context of roughening. 
Leaving the matter of preroughening aside, we will try to 
elucidate the effects of long-range elastic interactions on 
the roughening process. 

The paper is organized as follows. First we will intro- 
duce elastic interaction between surface defects and sug- 
gest some simplifications to make the problem tractable. 
Then we present the details of our discrete solid-on-solid 
model allowing for long-range step interaction. We will 
show the results of our extensive Monte Carlo simulations 
and interpret the effects. 



II. STEP INTERACTION 

Elastic step interaction on the surface of a semi-infinite 
crystal can be described i n || tc j| 'ni > ] of clastic force dipoles 
located at the step edgesOllalllUlS. Knowing the Green 
function Gy for an infinite elastic half-space one is able 
to calculate the elastic displacement field Ui{r) from a 



given force density fi{r) 

u,ir)= f dVG., (r-r')/, (r'). (1) 



The elastic energy Ec\ becomes 

E^, ^ ~ I <fru,{r)f,{r) 

d\d\'G,,{r-r')fj{r')U{r). (2) 

Using that forces fi{r) are present only in the vicinity 
of a step and that the mono-pole moment at the step 
vanishes, we can rewrite the energy using force dipole 
densities qik{r) as the next term in a multi-pole expan- 
sion 



Eel 



d?rd\' q,k{r') qu{r)dkdiG.,j{r - r'). (3) 



Using symmetry arguments one can determine two 
types |©f force dipoles that are considered to be present at 
a steplla. One type involves in-planc forces perpendicular 
to the step, the other arises from forces orthogonal to the 
crystal surface. Due to the structure of the Green func- 
tion, dipole tensors involving forces orthogonal to the sur- 
face show a behavior different from those involving only 
in-plane forcealjilj. The former lead to attractive or re- 
pulsive interaction depending on the signs of the steps, 
the latter produce a sign independent behavior, which 
is strictly repulsive. There are materialalj'El where the 
sign dependent contributions are small compared to the 
step repulsion caused by in-plane forces, and we will re- 
strict our model to the case, where we can neglect sign 
dependence of the steps. 

In the case of isotropic linear elasticity the half-space 
elastic. Green function Gij (r) can be written in a simple 
formlla 



i{(l-.)5,+a^} (4) 



■kE r 

where i and j are restricted to in-plane coordinates. 



2 



For a step stretching in y-direction one would assume 
the force dipole tensor at the step to be of the type ~ 
SixSjx- This means that the interaction between two hne 
elements will depend on their orientation. 

In the case of two parallel steps, a distance d in y- 
dircction apart, the interaction energy density w (per 
area squared) can be computed by evaluating the inte- 
grand from eq. (||) for two interacting force dipoles of the 
type qij = 5.,x5jx- It is given by 



w{r, ip) = ^ 



3 cos — 1 



a 2 + 15 cos** Lp — Ih cos^ ip 
1 — cr 



(5) 



where Lp denotes the angle between the radius vector r 
and the orientation of the dipole forces, which is given by 
ip = arctan(Ay/(i), and Ay is the distance between the 
dipoles in y-direction. The factor 7 is given by 



7 = 



l_oj 



(6) 



where Q = Q/a is the dipole moment (per unit step 
length) and Q would be the dipole moment assigned to 
a single atom at the step edge. 

Integrating the energy density for a configuration with 
two parallel steps at distance d, we state that the energy 
per unit length of the line is just 



W ■ 



. 1 „ 1 1 



where the interaction was limited to distances greater 
£. Note that the second term, which contributes to line 
energy, is negative for all possible Poisson ratios — 1 < 
a < 1/2. 

In order to make another simplification of the step- 
step interaction we compare the above result to the case 
of a scalar w ^ 1 jr^ interaction associated with isotropic 
dipoles ^ Sij, 



Ws. 



:alar 



1 1 

4735 + 2771, 



(8) 



c?2 

from which we conclude that the only difference in this 
specific geometry is a change in the line energy, which 
is mainly due to contributions from short range interac- 
tions. 

Because we aim at showing the effect of long-range 
interactions on the thermal roughening process, we ne- 
glect the angular dependence completely and assume that 
the dipole moments are isotropic. This leads to a simple 
isotropic 1 / r'^-interaction between force dipoles. Further- 
more this ensures that the elastic contribution to the step 
energy is positive. 



III. MODEL DESCRIPTION 

Within the framework of a solid-on-solid model we de- 
scribe the crystal surface by a simple height field of inte- 
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FIG. 1: The numbers of dipole charges assigned to a lattice 
sites is proportional to the accumulated absolute height dif- 
ference corresponding to eq. 



ger numbers h. Like in a common SOS model, overhangs 
are forbidden. Instead of the usual surface energy term 
(summation over nearest neighbors) 



E 



\hi - hj 



(9) 



with coupling constant J and a = 1, 2 for the ASOS- 
model and the DGSOS-model respectively, we define an 
elastic step interaction by introducing a field of elastic 
dipole charges q. To every lattice site a dipole charge 
Qk proportional to the number of height differences to 
the four neighboring sites is assigned, i. e. site k carries a 
number of 



Qk 



(10) 



<»,i> 



charges. Figure |l| gives an example how charges are as- 
signed to a simple height field configuration. The elastic 
dipole charges interact, in consequence of eq. (^, via a 
modified interaction potential ^'(r), 



*(r) 



r.-3 
1 



if r > 1 
if r 



(11) 



where r is the in-plane distance between two lattice sites 
measured in units of the lattice constant. This gives rise 
to the elastic energy 



Wcl \ ^ -r , N 



(12) 



^,3 



where Vij is the distance between lattice sites i and j 
and Wcl can be adjusted to give the desired interaction 
strength. Note that the case i = j is not excluded from 
the summation. During our simulation this constant had 
the value Wd = 0.31, which gives a line energy of about 
0.5 per unit length for infinite range interaction. 



Later we also limit the range of interaction. For this 
purpose we introduce a cutoff-potential 5"; with cutoff 
length I 

Mr)-[ y if,;; (13) 

which vanishes for distances greater than I. 

For two straight steps of length L with distance d this 
elastic energy contribution consists of the self energies of 
the steps and the expected ~ d^^ step interaction term 

Eint-^. (14) 

The self energy contribution of a straight step can be 
adjusted to be the same as the line energy of a DGSOS 
model. Because in this model both the line energy and 
the step interaction originate from eq. (p^, their rela- 
tive amplitude is fixed and we can concentrate on the 
crossover from a local to a long-range model depending 
on the cutoff length I. Note, however, that other relative 
amplitudes can in principle be obtained by a different 
choice of 5'(r) at small distances r < 1 in eq. ([ll]). 

The simulation is carried out on a square lattice of size 
i2 ^ 64 X 64 to 128 x 128. In order to calculate the differ- 
ence in energy for every metropolis Mnnte Carlo trial, we 
apply a multi-grid scheme based ono which has akeady 
been applied successfully to submonolayer epitaxyE3. 

This cuts down the computational costs from order L'^ 
to order L^\og{L) for each time-step, which has to be 
multiplied by an additional factor of L^, for the number 
of time-steps the system needs to equilibrate. Without 
the use of the multi-grid scheme the computational costs 
would not have permitted system sizes beyond L = 25. 
Still the system size, L < 128, is rather restricted and we 
are aware that the results should be accounted as qualita- 
tive rather than quantitative. However, computations on 
the DGSOS and ASOS models at L = 128, which we did 
for comparison, give transition temperatures Tr w 1.5 J 
and Tr » 1.25 J respectively, which agree reasonably well 
with known resultscS. 



IV. RESULTS 

A. Height Correlation Function 

We determine the roughening temperature Tr from the 
behavior of the height-height correlation function. Below 
roughening, T < Tr, the interface is macroscopically flat, 
i. e. the height-height correlation function 

G(r) = ([MO)-M^)]') (15) 

approaches a finite value in the limit r ^ oo. To be 
more precise, the correlation length ^ is finite and the 
interface has a characteristic width. Approaching the 
transition temperature the correlation length increases 
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FIG. 2: Height-height correlation function without cutoff. 
The correlation function saturates for small temperatures and 
shows logarithmic behavior for T > Tr,. The first straight line 
gives an estimate of Tr ~ 2.8. 



and diverges at T = X«. For T > Tr the correlation 
function G(j) diverge^3 according to the conventional 
theory of the roughening transition, 

G{r) ^ K{T)\o'gr, (16) 

with an amplitude K{T) depending on the temperature. 
Plotting G{r) vs. logr one could determine at what tem- 
perature the correlation length ^ diverges and the graphs 
approach a straight line. 

In a finite system with periodic boundary conditions, 
however, the correlation length ^ cannot exceed the sys- 
tem size L, the height-height correlation function G{r) 
saturates for T > Tr as well. In order to overcome this 
finite-size effect, we will use an approach similar to the 
one used inEj. In order to keep the argument simple we 
only consider correlations along the main directions of 
the lattice and replace r by x. 

As the limiting behavior of G{x) for periodic boundary 
conditions has to be a periodic function that behaves like 
the logarithm for distances <€i L, we define a "periodic 
logarithm" by means of Fourier analysis. In order to 
avoid the singularity at r — > we start with 

v{x) = max(log(x), 1) (17) 



4 



CUTOFF LENGTH = 1 



64x64 




FIG. 3: Height-height correlation function with a cutoff 
length 1, i.e. only charges on nearest-neighbor-sites interact. 
The first straight line gives an estimate of Tr ~ 0.9. 



and the integral-Fourier or, using symmetry arguments, 
the cosine transform 

oo 

v(k) = — / cos(kx)v(x)dx. (18) 

71- J 



Making use of these Fourier components we define the 
L-periodic function 

= cos — x^^, (19) 

11=1 V / V / L 

which is a discrete back transform averaged over unit 
distances. For convenience we define 

V{x)^Vl{x)~Vl{L/2) (20) 

and plot G{x) vs. V{x). Fig. || shows the correlation func- 
tion for the case of the full l/r^ interaction. At a temper- 
ature of about T = 2.8 the graph becomes straight, indi- 
cating the roughening transition. Restricting the elastic 
dipole charge interaction to distances < 1, the graph of 
the correlation function becomes straight at a lower tem- 
perature T = 0.9, see Fig. ^. 

From the Kosterlitz-Thouless theory of the roughen- 
ing transition, the slope of the correlation function is 
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FIG. 4: Slope K(T) vs. temperature. Estimation of Tr using 
the universal value from conventional roughening theory gives 
Tr = 2.8 for infinitely ranged interaction and Tr = 0.9 for 
cutoff length 1. 
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FIG. 5: Roughening temperature Tr vs. inverse cutoff length 
l/l. Even at / = 3 the roughening temperature Tr « 1.8 is 
well below the infinite range potential value. 



expected to assume the universal value K{Ty{) — 2/tt'^. 
Plotting slope vs. temperature one can also obtain an es- 
timate of the roughening temperature, see Fig. ^. From 
this we obtain identical estimates for the two cases with 
or without cutoff. 

From this we conclude that the system with long-range 
interaction has a much higher transition temperature 
compared to the model with interaction cutoff, the rough- 
ening temperature changes by a factor ~ 3. 

Note that the roughening temperature increases grad- 
ually with the cutoff length, see Fig. ^j. Even at / = 3 the 
roughening temperature Tr « 1.8 is still well below the 
value for infinite range interaction. The increase of the 
roughening temperature is not a next nearest neighbor 
effect. 



B. Energetic Scales 

One might argue that increasing the range of the inter- 
action potential just changes the relevant energetic scale. 
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FIG. 6: Average energy per unit surface E vs. temperature 
T. Average energy for the cutoff-potential is strictly higher 
in comparison to the long-range case. 
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FIG. 7: Scaled average energy per unit surface E vs. scaled 
temperature T/Tb.- For T/Tb. > 1 the scaled data collapses 
onto a single graph. 



However, the energetic scales one is usually tempted to 
think of, i.e. the energy of a straight step or single kinks 
on such a step, do not change by more than 36%. The 
straight line energy per unit step length increases from 
wi = 0.39 to 0.53, the corresponding kink energy changes 
from ~ 0.88 to 1.18. In the low temperature regime, 
the energy of one single ad-atom on a flat crystal sur- 
face is the important energetic scale, which changes from 
Wa = 2.53 to 2.69, an increase by no more than 6%. 

It should be noted that the main contribution to the 
change of these energetic scales comes form short-range 
interactions. Using a cutoff length of Z = 3, the straight 
line and kink energies are only about 6 — 7% below the 
the full potential value, whereas the single ad-atom defect 
energy deviates by no more than 0.05%. 

From the change of these energetic scales one usually 
would expect an equal increase of the roughening tem- 
perature. One would hesitate, however, to make these 
changes responsible for an increase of the roughening 
temperature by a factor of ~ 3. 



C. Average Energy 

Comparing the average energy E of the system com- 
puted with and without restriction of the charge interac- 
tion range, one clearly sees that the energy for the non 
restricted interaction always stays well below the graph 
of the restricted system, see Fig. ^. For high tempera- 
tures the average energy E goes linear with temperature 
T, indicating that the heat capacity becomes constant. 

The range of the interaction potential only affects the 
behavior below the transition temperature. Above the 
roughening transition all details of the interaction are 
combined into one single parameter, the roughening tem- 
perature Tr. Accordingly the scaled graphs E/T^ vs. 
T/Tr coincide for T/Tr > 1, see Fig. ^ 

The decrease in average energy of the system using 
long-ranged interaction coincides with a smaller number 
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FIG. 8: Extra surface, i. e. number of broken bonds, vs. tem- 
perature T. Fewer defects are created when no cutoff is used. 



of broken bonds, see Fig. ||. The number of deviations 
from a facet or the step length is smaller compared to 
the system with interaction potential cutoff. 



D. Defect Correlations 

Restricting the surface height to { — 1,0,-1-1}, one may 
talk about a defect wherever the height deviates from 
the average height 0. Then one can analyze correlation 
between these defects, i.e. the thermal average of 



{[h^r')-h^{r' + r)]^)r 
{h'{r'))l, 



(21) 



Scaled like this, the defect correlation will approach the 
value 1 for large distances r. At low temperatures repul- 
sion between the defects causes the graph to fall below 
value 1 at midrange distances and ends well above value 
1 at distance r = 1, because contact between equal de- 
fects is favored due to what might be called surface or 



6 




YL \ \ \ 

1 2 4 8 16 32 

X 



FIG. 9: Defect correlation at T = 0.5 for both long-range 
interaction and cutoff. The long-range interaction causes a 
stronger repulsion gap (here near x = 2), which means that 
the defects prefer to be separated. This anti-correlation effect 
is responsible for a strong decrease in entropy. 



line energy. Increasing the temperature, this repulsion 
gap will become smaller and vanish eventually. 

Fig. ^ shows the defect correlation for both long-range 
interaction and cutoff for identical temperature. Whereas 
for long-range interaction the gap is still present, it has 
already vanished from the system with cutoff. The pro- 
nounced repulsion gap in the case of the infinite range 
interaction means that a single defect or island avoids 
being close to other defects. This cuts down the number 
of favorable configurations and thus reduces the entropy 
contribution to the free energy for given density of defects 
n. 

For the following argument we will assume that the 
main result is a reduction of entropy by some factor 



a < 1 , whereas the average energy at given n remains un- 
changed. In a rather simplified picture we can then write 
the free energy like = E{n) — TaS{n), where n de- 
pends on temperature T and is determined by dF/dn = 
0. In this picture the free energy Fa{T) of the system 
with reduced entropy at temperature T has the same 
properties as the original system at a lower temperature 
aT. Thus if the original system had a roughening tem- 
perature Tr the transition temperature of the system 
with reduced entropy will increase to Tr = Tpj/a. 



V. CONCLUSION 

In summary, we have presented a model which con- 
tains the essential effects of long-range elastic repul- 
sion between steps on a crystal surface. We conclude 
that correlations due to these long-range interactions can 
strongly increase the roughening temperature in solid-on- 
solid models, mainly by a reduction of the entropy. Since 
defects prefer to exist in secluded areas, the number of fa- 
vorable configurations and consequently the entropy con- 
tribution to the free energy is diminished, leading to an 
increase of the roughening temperature. Our simulations 
suggest that the type of transition remains the same, al- 
though a rigorous proof lies beyond the scope of this type 
of Monte Carlo approach. 
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